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A  NUMERICAL  METHOD  FOR  THE 
TIME- DEPENDENT  TRANSPORT  EQUATION 

by 

R.  D.  Richtmyer 
ABSTRACT 

A  finite  difference  method  has  been  devised  for  the 
numerical  solution  of  the  Integro-dlff erentlal  equation 
of  neutron  transport  in  a  sphere.   The  method  is  applied 
here  to  one-group  problems  with  isotropic  scattering,  al- 
though some  degree  of  generalization  is  possible.   Special 
features  of  the  method  Include  integration  of  the  equations 
along  the  neutron  trajectories  in  space- time  and  special 
treatment  of  exceptional  net  points  which  lie  near  the 
boundary  of  the  relevant  region  of  phase  space.   The  method 
is  the  outcome  of  a  series  of  Investigations  initiated  by 
a  suggestion  of  Von  Neumann  in  19^8. 

Numerical  tests  of  the  method,  made  with  the  Unlvac 
at  New  York  University,  are  described.   They  deal  with  the 
problem  of  a  homogeneous  bare  sphere  subject  to  a  singular 
initial  condition.   In  connection  with   those  tests,  codes 
were  divised  to  provide  accurate  numerical  Integration  of 
the  Integral  equation  for  the  normal  mode  functions. 
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A  NUMERICAL  METHOD  FOR  THE  TIME- DEPENDENT 
TRANSPORT  EQUATION 

Introduction 


A  finite  difference  method  for  solving  the  time- 
dependent  transport  equation  In  spherical  systems  was 
devised  by  John  Von  Neumann,  Herman  H.  Goldstlne  and  the 
writer  (unpublished)  In  1948  In  connection  with  work  of 
the  Los  Alamos  Scientific  Laboratory  and  was  used  in  ex- 
tensive machine  calculations  in  19^9  and  1950.   The  essen- 
tial feature  of  the  method  was  the  use  of  the  quasi-cartesian 
coordinates  described  below.   In  1952  and  1953  the  method 
was  further  improved  and  simplified  by  integrating  along 
the  neutron  trajectories  In  space- time;  preliminary  tests 
of  the  method  (unpublished)  in  this  form  were  made  on  the 
Los  Alamos  computer  ("Maniac").   During  the  past  year  and 
a  half  further  tests  have  been  made  on  the  Univac  at  New 
York  University  for  a  simplified  problem  and  during  the 
course  of  these  tests  some  further  substantial  Improvements 
of  the  method  have  been  devised,  to  provide  more  accurate 
treatment  of  exceptional  net  points  near  the  boundary  of 
the  system. 

The  latest  version  of  the  method  and  the  Univac  tests 
will  be  described  here. 
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For  practical  problems  the  time-dependent  version  of 
Carlson's   S   method  has  been  generally  preferred  in 
recent  work,  because  of  its  smaller  demands  on  machine 
speed  and  capacity;  but  the  method  presented  in  the 
present  report  is  of  Interest  because  of  its  directness 
and  simplicity. 
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II.   The  Initial  Value  Problem 


To  simplify  discussion,  attention  will  be  confined 
to  the  simplest  possible  version  of  the  problem  by  making 
the  following  assumptions: 

1.  All  neutrons  have  the  same  speed  ("one-group 
treatment") . 

2.  Scattering  is  isotropic. 

j5.   Cross-sections  and  other  material  parameters 
are  time- independent. 
It  will  probably  be  obvious  to  the  reader  that  some  degree 
of  generalization  is  possible  and  necessary  for  practical 
problems. 

Let  r  =  radial  coordinate;  m-  =  cosine  of  the  angle 
between  the  radius  vector  and  the  velocity  vector;  t  =  time; 
v  =  neutron  speed;  6  =  6(r)  =  collision  probability  per 
unit  path  length  (also  called  the  macroscopic  cross-section 
or  inverse  mean  free  path);  1  +  f  =  1  +  f(r)  =  average 
number  of  neutrons  emerging  from  a  collision,  by  scat- 
tering and  fission;  il/{v,\i,t)    =   number  of  neutrons  at 
time  t  per  unit  volume  of  space  at  position  r  per 
unit  solid  angle  about  direction  |i.   Then  the  transport 
equation  is 

=  6  ^  J     dn'TA(r,n',t). 
^   -1 
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This  equations  Is  to  be  satisfied  for  0<r<a,  -1<  m.<1, 
0  <  t,  where  a  Is  the  (outside)  radius  of  the  system;  the 
equation  is  to  be  taken  together  with  an  initial  condition: 

(2)         ^(r,|j.,0)  =  •^^(r,!^)   (given), 
and  a  boundary  condition: 

(5)         T^(a,pL,t)  =0  for  -1  <  n  <  0. 

If  the  reader  so  desires,  he  may  regard  the  requirement 
that  all  terms  of  (l)  be  bounded  as  r  ->  0  as  an  addi- 
tional boundary  condition  holding  at   r  =  0. 

Many  properties  of  this  initial  value  problem  are 
well-known,  at  least  if  the  functions  C5(r)  and  f(r) 
are  sufficiently  well-behaved  --  for  example  piecewise 
continuous  .   Some  of  these  properties  are: 

1.  A  unique  solution  exists  for   t  >  0   if  the 
initial  function  V  (r,|jL)  is  piecewise  con- 
tinuous . 

2.  The  solution  depends  continuously  on  the 
initial  function  in  the  sense  of  the  maximum 
norm. 

^.   The  solution  cannot  grow  faster  than  at  a 
certain  exponential  rate,  naunely: 

(4)   U(r,^i,t)|  <  e^  Max   i^^(r,^;)| 

o<r<a   ^ 
-1<H<1 


•Piecewise  continuity  is  taken  to  Include  possession  of  one- 
sided limits  at  each  discontinuity,  hence  boundedness. 
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where 


B  =  v[Max  6(r)]  [Max  f(r)] 
o<r<a       o<r<a 

4.  The  solution  Is  continuous  In  r  and  p.  for 
t  >  2a/v  even  If  the  Initial  function  Is  dis- 
continuous. 

5.  The  general  solution  of  the  Initial  value 
problem  can  be  expressed  In  terras  of  normal 
mode  solutions  and  a  residuum  as  follows: 

(5)    T^(r,n,t)  =  X:(n)  ^n   >^^"^^>^^)e°'  ""  ^   +   d(r,n,t) 

where,  for  each  n  =  0,1,2,...,   a^'^'  and  i/   ^"' 

are  a  constant  and  a  function  characterizing  the 

n   normal  mode;  the  constants   A   and  the 

n 

residuum  d(r,p.,t)   depend  on  the  Initial  func- 
tion but  In  any  case  can  be  so  chosen  that  as 
t  — ^>  00  ,   d(r,|j.,t)   decreases  faster  than  any 
normal  mode. 
(This  last  Is  a  conjecture  based  on  the  studies  of 
homogeneous  bare  slabs  and  spheres  by  Lehner  and  Wing  -- 
see  article  by  Lehner  In  Communications  on  Pure  and 
Applied  Mathematics,  vol.  IX,  1956.   For  the  slab  the 
s\immation  in  (5)  is  a  finite  one  and  the  residual  term 
is  generally  necessary.   For  spherical  systems  the  summa- 
tion is  probably  an  infinite  one  and  the  residual  term 
may  not  be  necessary.) 

Clearly  this  Initial  value  problem  is  suitable  for 

attempts  at  numerical  solution. 

-  8  - 


NYO-7696 

Physics 


III.   The  Quasl-Carteslan  Coordinates 

In  place  of  r,  ly,      we  take  as  independent  variables 
the  quantities 


X  =  rp,  =  r  cos  "^  , 


(6) 


y  =  r  /l-p.    =  r  sin  \T , 


where  we  have  written  |j.  =  cos  '^  . 

These  are  not  cartesian  coordinates  In  the  ordinary 
sense  because  r  and  \^    are  not  polar  coordinates  In  the 
ordinary  sense.  '^  Is  the  angle  between  the  radius  vector 
and  a  fixed  direction.   In  terras  of  the  straight  line  on 
which  a  neutron  is  moving  at  time  t,   y  is  the  shortest 
distance  from  this  line  to  the  origin  and  x  is  the  dis- 
tance along  this  line  measured  from  the  point  nearest  the 
origin  to  the  present  position  of  the  neutron,  the  posi- 
tive direction  being  taken  as  that  of  the  neutron's 
motion.   Thus,   x  and  y  are  cartesian  coordinates  in  a 
plane  which  shifts  its  orientation  abruptly  when  there 
is  a  collision. 

With  these  variables,  and  calling  ^(x,y,t)  =  V(r,cos i^t] 
the  transport  equation  becomes 

C^)   (^3^  +  jl  ^  ^)  *(x,y,t)  =  6  ^  (D  {/F^,t), 
where  $(r,t)  is  the  spatial  density  given  by 

IT 

(8)        <I)(r,t)    =  27r  /     ¥(r  cos  if,    r  sin  ir',t)sin  x/d  ^^ 

o 

or 
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(9) 


$ 


(r,t)  =^     /  dx  f(x,  /r^-x^  ,t) 


-r 


(7)  is  to  be  satisfied  for  x^  +  y^  <  a^  ,  y  >  0,  t  >  0; 
the  Initial  condition  is 


*(x,y,o)  given,   x^  +  y^  <  a^  ,  y  >  C  ; 


the  boundary  condition  is 


(10)   ^(x,y,t)  =0  for  x"^  +  y^  =  a^,  x  <  0,  y  >  0 
or 


^(-  /a^-y^  ,  y,t)  =0  for  0  <  y  <  a. 


There  is  no  boundary  condition  at  the  origin. 

The  transport  equation  holds  in  the  shaded  semi- 
circular region  of  figure  1  and  the  path  of  Integration  in 


Figure  1 

(8)  or  (9)  is  a  semicircular  arc  passing  through  the  point 
x,y.   In  this  figure  the  point  representing  a  neutron  moves 
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horizontally  to  the  right  with  speed  v  until  a  colli- 
sion occurs,  at  which  time  It  jumps  to  some  other  point 
on  the  same  semicircle  (say  x'y")  from  which  it  then 
continues  to  move  horizontally  with  speed  v.   If  a 
neutron  escapes  from  the  system  It  does  so  by  crossing 
the  right-hand  quarter  circle  of  the  boundary.   The 
boundary  condition  of  no  Incident  flux  is  that   ^  =  0 
on  the  left-hand  quarter  circle  of  the  boundai^. 
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IV.   The  Difference  Equations 

A  rectangular  net  of  points  In  the  x,y,t  space  Is 
chosen  by  defining 

x^  =  lAx        1  =  -1,-1+1,- •• ,0,1,- •', I 
yj  =  JAy        J  =  0,1,- ••,J 
t^  =  nAt        n  =  0,1,- • • 

where  the  Intervals  are  supposed  so  chosen  that 
vAt  =  Ax  and   lAx  =  JAy  =  a.   f(xj^,y  ,t")   will  be 
denoted  by  *J  j  •   <D(r,t^)   will  be  denoted  by  <l)'^(r). 
Furthermore  an  Interval  Ar  Is  chosen  for  tabulation  of 
$  (r)   and  we  will  denote  <I)"(pAr)   by  $".   We  also 

Xj^+y   by  r^  ..   Fractional  indices  will  also 
be  used  where  convenient. 

Because  of  the  relation  vAt  =  Ax,  the  space-time 
point  (^±+i'y  i'^^      )   lies  on  the  same  trajectory  as 
the  point   (x. ,y,,t  ),  and  we  can  approximate  the  deriva- 
tive terms  by  writing 

^^i+l/2'yj'^      ' 
and  the  transport  equation  is  approximated  as 


•The  trajectories  or  characteristics  are  the  lines 
x  -  vt  =  constant,  y  =  constant. 
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(11) 


Ax  2 


^I+f  ^n+x/2  / 


where  ^^'^  '     (r)   is  an  approximate  value  of  (i)(r,t)  for 
time  (n-t-p)At,  to  be  obtained  by  an  Iterative  procedure 
described  below.   If  the  ^'^  '         are  available,  equation 
(11)  gives  the  unknowns   f     explicitly  at  all  net  points, 
except  that  if   (l+l,j)   denotes  the  left-most  point  on  the 
j   horizontal  line,  inside  the  semicircle,  like  the  point 
B-,  In  figure  2,  the  quantity   f"  .   needed  in  (11 )  is  not 
available.   But  the  boundary  condition  says  that  1.   vanishes 
at  point   P-,   of  space-time  at  which  the  trajectory  through 
(i+l,J,n+l)  enters  the  system,  at  radius  a.   In  this  case 
we  replace  Z^x  in  (ll)  by  /a  -y  ^.  -  x,  ,  ,  which  is  the 
distance  from  the  point   (i+l,j)   to  the  boundary  along  the 
line  y  =  y.;  we  also  replace   f^  .   by  zero;  then  (ll) 
gives   ^      as  for  other  net  points  in  the  system. 

To  explain   more  fully  the  treatment  of  exceptional 
points,  we   refer  to  figure  2,  which  is  a  perspective  drawing 
showing  certain  features  of  the  three-dimensional  net.   The 
line  segments   A, Bp,  ApB,,  etc.  are  trajectories  or  character- 
istics connecting  net  points  at  time   t   with  those  at 
time   t'^'^''".   The  quantity  ^'^'^''^    ^^-1/2  1^   ^'^  equation 
(11)  represents  generally  an  average  value  of  <t>(r,t)   along 
such  a  segment  and  is  taken  as  the  average  value  of   $  at  the 
two  ends  of  the  segment.   For  a  segment  like   A,B2   this 
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involves  a  value  of  <^^{v)    at  one  end  and  of  $""*'■'•  at  the  other. 


planes   t  =  t   and  t  =  t 


Figure  2 
The  special  short  segments  like  P-iB-,   and  PpP^ 
require  a  modification  of  the  procedure.   These  segments 
terminate  on  the  plane  t  =  t     but  originate  at  points, 
namely  P.,   and  Pp  respectively,  which  lie  between  the 

We  first  obtain  a  value  of 
$  at  the  early  end  of  such  a  segment  by  linear  interpolation 
in  time;  then  the  quantity  o"^  '    in  equation  (ll)  is  taken 
as  the  average  of  the  values  of  $  at  the  two  ends  of  the 
segment.   For  example,  in  the  equation  to  determine  Y 
at  point  P-,  (see  figure  2),  the  right  hand  side  of  (11 ) 
is  taken  as 


C5^  [|  ^(Pg)  +  I  <^iP^)]       . 


where  <I>(Pp)   has  been  obtained  in  turn  by  linear  inter- 
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polatlc*^  b«=tween  <j>{A  )      and  <i)(Bg).   The  net  result  is 
that  (^"■'■■'■/^(r,   /„  ,)  In  equation  (ll)  has  been  Inter- 
preted  as  a  value  obtained  by  two-way  interpolation  (with 
respect  to  x  and  t)  among  the  values  of  $  at  points 
A  ,Bq  and  P^   to  give  an  approximation  to  $  at  the 
midpoint  of  the  segment  F^P-. 

A  further  refinement  of  this  equation  (the  equation 
to  determine  y""*""^  at  point  P^)  is  that   y"  .   in  the  left- 
hand  member  is  similarly  replaced  by  "ii^o)      which  is  a 
value  obtained  by  linear  interpolation  between  f(A,.)   and 
^(B  )  --  that  is,  between   '^^  .   and   ^^"t"'-  .   (This 
requires  that  the  equations  be  solved  in  order  of  increasing 
i,  so  that  Y?t   will  be  available  when  we  are  computing 
^(P.)).   Lastly   Ax   is  replaced  by  the  distance  between 
the  points   B   and  P   --  i.e.  by  /a  -y.  -  x. . 

When  the  ^.  ,   are  Known  at  all  net  points,  the 
$^"''   are  computed.   Equation  (8)  or  (9)  shows  that  $" 
is  given  by  a  line  integral  of  ^'^'^       around  a  semicircle 
of  radius   pAr   in  the  x,y  plane.   Generally  this  semi- 
circle will  pass  through  few,  if  any,  of  the  net  points, 
and  interpolation  is  required.   The  machine  calculates 
<l>_   by  applying  the  trapezoid  rule  to  equation  (9),    using 
2"  Ax  as  the  increment  of  x  for  this  purpose  and  obtaining 
the  values  of  the  integrand  at  the  required  abscissas   x 
by  linear  interpolation  with  respect  to  x  and  y. 
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The  Iterative  procedure  for  solving  (ll)  is  as  fol- 
lows:  At  the  beginning  of  a  complete  cycle  of  the  calcula- 
tion, the  ^^    .     are  known  and  a  table  of  ^'^  is  available. 

1  >  J  P 

The  quantity  ^^'^■^^^    ^^1+1/2  1^   °"  ^^^   right  of  (ll)  is 

first  approximated  by  simply  4>  (r.  ■>  /^  .) ,    obtained  by 

linear  interpolation  in  the  table  of  *   at  r  =  r,  .,  /„  .. 

P  ±+1/ d , J 

Equation  (ll)  is  then  solved  to  give  preliminary  values 

of  ^7"'"..   Prom  these,  provisional  values  of  $""*"  are 
i,J  P 

obtained  by  the  trapezold-rule  integration  described  in 

the  preceding  paragraph.   A  second  calculation  (iteration) 

is  then  performed,  in  which  i>^'^    '    (r)   is  taken  as  the 

mean  of  <t>"(r)   and  (^^    (r).   This  iteration  results  in 

improved  values  of  \,    ,    ,  from  which  improved  values  of 

^        are  obtained. 
P 

Experience  indicates  that  generally  further  itera- 
tions are  not  worthwhile  --  using  a  finer  net  pays  off 
more  rapidly.   This  is  presumably  because  the  first 
iteration  decreases  the  truncation  error  formally  from 
(?^(At)   to  (5'((At)  ),  and  further  iterations  change  only 
the  coefficient  of  the  error  term  but  not  the  power  of 
At. 
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V.   The  Trial  Calculation 


The  method  was  applied  to  a  simplified  problem  of 
conceivable  astrophysical  interest  but  chosen  mainly  to 
illustrate  the  method.   It  applies  to  photons  rather 
than  neutrons;  hence  the  speed  is  called   c   instead  of 
V.   A  short  burst  of  light  (as  from  a  variable  star)  is 
emitted  at  time  t  =  0  from  a  point  source  situated 
at  the  center  of  a  large  homogeneous  spherical  cloud 
of  purely  scattering  material  (f  =  0)  which  is  assumed 
to  scatter  isotropically  and  without  polarization 
effects.   Required  is  a  curve  of  intensity  vs.  time 
for  the  light  emitted  from  the  cloud. 

Several  calculations  were  made  for  a  sphere  of  2 
mean  free  paths  radius  and  one  calculation  for  a  sphere 
of  4  mean  free  paths  radius.   Some  of  the  results  are 
given  graphically  in  figures  3A,  3B,  4,  5,  6.   Choosing 
the  mean  free  path  and  mean  free  time  as  units  of  length 
and  time,  we  have: 

Case  I      6=1,   v=c=l,   a=2. 
Case  II     6  =  1,   v=c=l,   a.  =  4. 

In  the  first  two  attemps  to  calculate  case  I  (I  =  J  =  10, 
162  interior  net  points;  and  I  =  J  =  15,  about  550  interior 
net  points),  the  problem  was  treated  in  a  straightforward 
way  as  an  initial- value  problem:   at  n  =  0  (t  =  O)  Y^, 
was  set  equal  to  zero  except  at  the  central  net  point 
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1  =  J  =  0,   where  it  was  arbitrarily  set  equal  to  0.1. 
The  photons  were  thus  initially  concentrated  into  the 
immediate  neighborhood  of  the  point  source  and  then 
allowed  to  move  in  accordance  with  the  transport  equa- 
tion.  This  treatment  was  unsuccessful,  except  in  a 
rough,  qualitative  way,  because  of  the  violently  dis- 
continuous nature  of  the  distribution  at  early  times. 
Spreading  out  the  initial  distribution  over  5  or  ^  net 
points  was  also  tried,  but  produced  little  improvement. 

The  procedure  was  then  modified  as  follows:  ^^ ,  was 
taken  identically  zero  for  n  =  0,  and  a  source  term  was 
introduced  into  the  transport  equation  during  the  interval 
0  <  t  <  a/c  while  the  initial  light-pulse  from  the  star 
was  travelling  outward  through  the  sphere,  the  source 
being  located  at  position  r  =  ct  at  time  t  and 
representing  photons  emerging  from  their  first  scatter- 
ing collision.   In  other  words,  the  first  collision  was 
treated  analytically  and  subsequent  ones  by  the  nvunerical 
calculation.   The  Initial  pulse  travels  out  as  an  expand- 
ing spherical  shell  with  a  number  of  photons  per  unit 
area  of  the  shell  proportional  to  e~  ^/r    .      To  allow 
for  a  shell  source  proportional  to  this,  it  is  necessary 
only  to  add  the  quantity   Ae'^^'^/p^   to  ^^'^      after  the 
cycle  for  which  n+1  =  p.   Before  the  first  cycle  (t)^  was 
set  equal  to  a  value  chosen  to  represent  approximately 
the  photons  injected  into  the  system  by  the  shell  source 
during  the  time-interval  {0,-^lst)    (the  other  <^^     were 

-18- 


NYO-7696 

Physics 

zero) .   The  constant   A  was  taken  proportional  to   I 
so  that  when  a  problem  was  rerun  with  a  finer  net  the 
total  niimber  of  photons  introduced  during  the  interval 
(0,a/c)   should  be  the  same.   The  number  of  photons 
Introduced  In  the  p   cycle  by  our  procedure  Is  pro- 
portional to 

5$  r^ArAt  =  6$  (pAr)^Ar^  = 
p  p  *    '    c 

p  c   ^I' 

where  6$  =  Ae~^  /p  .   The  total  number  Introduced  is 
proportional  to 

o 

With  this  modification  the  functions  i>     and  *  which 
one  Is  trying  to  compute  still  have  discontinuities  during 
the  time  interval   (0,a/c),  but  at  least  they  are  bounded 
cind  do  not  require  Dirac  delta-functions  for  their  repre- 
sentation. 

The  anergent  flux  F  at  r  =  a  was  obtained  by 
trapezoid-rule  evaluation  of  the  integral  in  the  equation: 

Tr/2 
F  =  2ir  J         sin^yd^r^Ca  cos  1^  a  sini^,t). 
o 

The  flux  P  is  shown  plotted  on  a  logarithmic  scale 

as  function  of   t  for  case  I  in  figure  5a  and  for 
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case  II  in  figure  3b.   In  both  cases   F  Is  practically- 
zero  (as  It  should  be)  until  the  Instant  t  =  a/c  when 
the  original  pulse  emerges.   Thereafter  the  cloud  shines 
with  an  intensity  which  varies  with  time  as  shown.   At 
first  the  curve  of   F  vs.  t   shows  a  transient  behavior, 
but  later  it  settles  down  to  a  very  nearly  exponential 
decline. 

The  spatial  density  $  =  «|)(r,t)   is  shown  as  function 


Figure  4 
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of  r  at  various  times   t  for  case  II  In  figure  4.   In 
the  three  earliest  curves  the  discontinuous  radiation 
front  at  r  =  ct  is  clearly  visible.   The  oscillations 
behind  the  front  are  presumably  spurious,  having  been 
somehow  Introduced  by  the  truncation  error  of  the  finite 
difference  equations.   After  the  instant  t  =  a/c  =  2, 
the  distribution  Is  quite  smooth  and  settles  down  to  a 
steady  shape. 

We  have  no  exact  information  with  which  to  compare 
these  results  but  one  can  get  an  idea  of  the  accuracy 
attained  in  various  indirect  ways.   In  the  first  place 
there  are  several  internal  checks  which  can  be  applied 
to  the  calculation.   For  example,  the  total  number  N 
of  photons  in  the  cloud  at  an  Instant  t  can  be  obtained 
either  as  the  space  integral  of  the  space  density: 

a  p 
N  =  N^  =  4Tr  /  r^  dr  $(r)   , 
o 

or  as  the  phase  space  integral  of  the  phase  space 
density: 

N  =  Ng  =  27r  //  dxdy  y  *(x,y)   , 

this  last  Integral  being  extended  over  the  semicircular 

2    2    2 
region:   x+y<a,y>0.    N,   and  Np  should  of 

course  be  equal.   After  each  cycle  the  machine  evaluated 

N-,   and  Np   (essentially  by  the  trapezoid  rule)  and 

printed  them  out.   During  the  Interval  0  <  t  <  a/c,  while 
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the  discontinuous  front  was  moving  outward,  there  was 
generally  a  quite  large  discrepancy  between  N-,   and  N^ 
but  later  the  agreement  was  better.   For  the  calculation 
shown  In  figure  3a,  the  difference  between  N-,   and  Np 
was  approximately: 

15%  at   t  =  I  a/c  =  1 

5  ^  at   t  =  a/c  =  2 

1. 6 ^( average)  for  a/c  <  t  <  2a/c 
(I.e.  for  2  <  t  <  4) 

Considering  the  small  number  of  net  points  used  (18  on 
the  radius),  this  discrepancy  Is  perhaps  not  surprising. 
Furthermore  there  is  some  evidence  that  even  this  dis- 
crepaj^cy  was  due  in  part  to  poor  scaling  (loss  of  significant 
figures)  rather  than  truncation  error,  for  the  portion 
2— <t<3—  of  the  calculation  was  performed  after  re- 
scaling,  and  the  average  discrepancy  between  N,   and  Np 
was  only  0.4  "/'  for  2  J  <  t  <  3  f  • 

As  a  second  internal  check,  the  quantity 

T  =  N  +  /  Fdt 

should  be  constant,  by  conservation  of  photons,  for  t  >  a/c; 
i.e.,  after  the  source  term  is  no  longer  present  in  the 
equation.   If  Np   is  used  for  N  in  the  above  equation, 
it  is  found  that  T  is  indeed  quite  constant;  the  average 
deviation  of  T  from  its  mean  value  for  a/c  <  t  /  2  a/c 
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was  0.11%  for  the  calculation  shovm  In  figure  3a.   This 
suggests.  Incidentally,  that   N^   (integration  over  phase 
space)  provides  a  better  estimate  of  N  that  does  N, 
(integration  over  space). 

A  second  type  of  check  on  the  accuracy  resulted 
from  varying  the  net  spacing.   The  calculation  was  also 
performed  with  a  coarser  net  than  used  for  figure  3a; 
I  =  12  Instead  of  I  =  18.   The  results  agreed  well  with 
the  ones  shown,  except  In  the  immediate  vicinity  of 
t  =  a/c,  where  the  sudden  onset  of  the  luminosity  was 
slightly  less  abrupt  than  for   I  =  18.   For  t  >  a/c, 
the  two  curves  have  the  ssune  shape  to  an  accuracy  of 
about  0.5%,  but  the  one  for  I  =  12  is  consistently 
about  2.5%  lower  than  that  for  I  =  18, 
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VI.   The  Asymptotic  Form 

Further  evidence  of  the  accuracy  of  the  method  comes 
from  the  asymptotic  behavior  of  the  solution  for  large 
values  of  t.   The  normal  mode  expansion  (5)  Indicates 
that  for  sufficiently  large  t  the  fundamental  mode 
(the  mode  with  the  highest  value  of  a)   should  predom- 
inate.  Therefore  the  asymptotic  form  is 


*(x,y,t)  =  Y(°^x,y)e°^  °  ^ 
*(r,t)  =  <t>(°^r)e°'^°^t 


If  these  expressions  are  substituted  into  the  transport 
equation,  the  resulting  equation  (which  then  has  the 
character  of  an  ordinary  differential  equation  for  fixed 
y  if  $^°^   is  regarded  as  known)  can  be  solved  for 
f'°^   in  terras  of  <I)'°  .   If  the  result  is  substituted 
into  the  equation  (8)  defining  $^°  ,  one  obtains  an 
integral  equation  for  <i)^°^(r),   which  can  be  put  into 
the  form 

a        a 


(12)     >.a>^^^(r)=i  /   E(((5+^)|s-r|)$^°^(3)ds 


where 


and 


^   =  (6(1  +  f))-^ 


OD   -y , 
E(x)  =   /  ^ 
X     "^ 
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cj)  (r)   is  In  fact  that  odd  solution  of  (12)  corresponding 
to  the  largest  eigenvalue  A.    (To  derive  (12)  one  must 
assume  that   6  Is  constant:   In  our  calculation,   f  was 
also  constant;  namely,  f  =  0) . 

The  numerical  method  of  solving  (12)  that  was  used 
will  be  described  In  the  next  section.   The  procedure  is 
more  or  less  standard  and  so  accurate  that  the  results 
may  be  regarded  as  exact  for  purposes  of  comparison  with 
those  described  In  the  preceding  section. 

Examination  of  the  numerical  solution  shown  In 
figure  3a  shows  that  for  2  —  <  t  <  3  —  the  decay  Is 
quite  accurately  exponential  and  the  value  of  the  decay 
constant  Is  a  =  -.490^1^,  from  the  slope  of  the  graph,  by 
least-squares  fit  to  the  last  few  points,  as  shown' In  the 
figure.   We  therefore  suppose  that  the  asymptotic  form  of 
the  solution  has  been  fairly  well  established  by  t  =  2  — 
and  we  compare  this  form  with  the  correct  one  as  given 
by  the  Integral  equation  method.   The  correct  value  of 
a  Is  a  =  -.4925,  from  which  our  value  differs  by  less 
than  1/2%.   In  figure  6  the  radial  dependence  of  <j>(r,t) 
for  t  =  2.5  —  (from  one  of  the  early  calculations)  is 
compared  with  the  fundamental  solution  <I)^°'(r)   of  (12), 
suitably  normalized.   The  discrepancy  is  within  l'% , 
except  for  the  point  on  the  surface,  r  »  2,  where  it  is 
where  It  is  3  7o  . 


-26- 


NYO-7696 

Physics 


Similarly,  a  least-squares  fit  gives  a  =  -.1515 
for  a  sphere  of  4  mean  free  paths  radius  (a  =  4),  as 
compared  with  the  correct  value  of  a  =  -.1529  from  the 
integral  equation. 
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VII.   Iterative  Solution  of  the  Integral  Equation 

In  writing  the  Integral  equation  (12)  we  have  fol- 
lowed the  procedure  that  Is  customary  In  this  work,  of 
defining  $  (r)   for  negative  r  as  -(i)°(-r)   to  simpli- 
fy the  form  of  the  result.   For  the  numerical  work  It  Is 
more  convenient  to  work  with  the  range  0  <  r  <  a. 
Recalling  that  f  =  0,  and  taking  units  such  that 
6  =  1,  v  =  1,  we  have 

a 

r 

o 
where 


(15>.     ?\ci)(r)  =  /  K^(r,s)  o(s)ds         0  <  r  <  a 


(14)     K^(r,s)  =  I  [E((l+a)|r-3|)  -  E( (l+a) | r+s | ) ] . 

The  superscript  "(o)"  has  also  been  dropped  from  <:)(r). 

For  any  real  a  >  -1,  K  (r,s)   Is  a  symmetric  kernal 
of  Hllbert-Schraldt  type  (non-degenerate)*.   It  can  also 
be  shown  to  be  positive  definite.   Therefore  the  eigen- 
values of  (i;5)  form  a  positive  sequence   A   >  ^i  >  ^p  '  '  ' 
tending  to  zero;  each  Is  of  finite  multiplicity;  and 
the  corresponding  eigenfunctlons  may  be  taken  as  an  ortho- 
normal  set  over  the  Interval   (0,a).   Furthermore,  since 


*  For  the  general  theory  of  such  Integral  equations  see 
Courant-Hllbert,  Methods  of  Mathematical  Physics,  Chapter 
III.   The  extensions  of  the  theory  necessitated  by  the 
singularity  of  the  kernel  K  (r,s)  at  r  =  s  are  described 
In  section  9  of  that  chapte?. 
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the  kernel   K  (r,s)   Is  positive  for  0<r<a,   0<s<a, 
the  highest  eigenvalue  >   is  of  multiplicity  1,  so 
X  >  A^ ,  and  the  corresponding  elgenfunction  $(r)   is 
of  one  sign  (and  may  be  taken  as  >  0)  throughout  the 
interval   (0,a).   The  elgenfunctions  corresponding  to 
^1*^2*'''  '^'^^^   t»^  denoted  by  X^(r),  X^{r) ,    etc. 

The  problem  is  to  find  a  value  of  a  such  that  the 
highest  eigenvalue  of  (13),  A  ,  is  exactly  =  6~      =  1, 
and  then  to  find  the  corresponding  elgenfunction   <I)(r). 
If  one  has  a  general  procedure  for  finding  the  largest 
eigenvalue  and  corresponding  elgenfunction  of  an  equation 
like  (13)  for  any  given  value  of  a,  then  it  is  only  neces- 
sary to  Include  a  trial-and-error  method  of  adjusting 
a  until   A   =  1. 

The  well-known  iterative  pr^/cedure  for  finding  the 
largest  eigenvalue  is  the  following:   Let  h(r)   be  an 
arbitrary  positive  continuous  function  for  0  <  r  <  a, 
and  define 

g^(r)  =  h(r) 

a 
gp(r)  =  /  dsK^(r,s)gp_^(s) 

p  =  1,2,3,  ••• 


a 
P 


{/  ds[gp(s)]2  ] 


P  =  0,1,2, 
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Clearly  a   Is  Just  the  reciprocal  of  the  normalizing 
factor  for  g  (r) .   According  to  the  general  theory  any 
function  that  can  be  written  as   /  K  (r,s)h(s)ds,  with 
some  continuous  h(s),  can  be  expanded  In  the  elgenfunc- 
tlons.   This  applies  to  g-,(r),  g^Cr),  etc.   Therefore, 

gl(r)  =  c^$(r)  +  c^X^(r)  +  c^XgCr)  +  •••, 


and 


vi  =  [  %V^  -  =1 V -3%'^  — ] 


1/2 


On  dividing  * ,  It  is  seen  that 


-^ =  <t>(r)  +  •••  , 

p+1 


where  the  dots  Indicate  terms  containing  factors 
\P     ^2p       >.P     >2p 

^h  •  C^J  ■■■■{€)    -{>)   ;•••--, 


since  Aq  >  Aj^  >  Ag  >  "  '  >   all  these  terms  vanish  In  the 

limit  as  p  _>oo.  That  Is, 

lira   ^p^^'  /  \ 

p^oo    a  '*'\  / 
P 


c   cannot  vanish,  because  h(r)  >  0,  therefore  gn(r)  >  0, 
therefore  c   =  /g,(r)$(r)dr  >  0  because  $(r)  >  0. 
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a 
Similarly,  11m  -^^  =   A^. 
p->oo    p 


VIII.   Improved  Iterative  Method 

There  are  many  well-known  schemes  for  accelerating 
the  convergence  of  this  Iterative  procedure,  and  the  one 
used  in  our  calculation   is  the  following.   Let   |j..  , (Xg  >  *  "  '  M-m 
be  N  real  constants,  whose  values  are  to  be  decided  on 
later.   Given  a  function  g  (r),  we  define 


gp+l^^^  =  ^^   K^(r,s)gp(s)ds  -  ^ilgp(r) 
(15)     gp^2(^)  =  /^  V^'^^Vl^"^^"  "  ^^2^1^^^ 


a 
Then  if   g  (r)   is  of  the  form   /   K  (r,s)h(s)   for  some 

P  o   ^^ 

h(s)   we  can  again  expand: 


gp(r)  =  c^O(r)  +  c^X-^(r)  +  ••• 

and  clearly 

Ep+^M    =  <^Qpi\)<^i^)    +  c^F{A^)X-^(r)  +  ••• 
where  P(>)   is  the  polynomial 

F(a)  =  (A-^L-^)(A-H2)•••(>^-^iN^■ 

« 
essentially  that  of  Flanders  and  Shortley,  J.  Appl.  Phys 

vol.  21,  p.  1326  (1950). 
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What  is  desired  is  to  choose  the  |j.,  ,  •  •  •  ,  |j.„  so  as  to 
make  the  ratios   P(>.j^  )/F(  A^) ,  Fi-K^)/F{-h^) ,    etc.,  all  as 
small  as  possible.   We  can  achieve  essentially  this  by 
use  of  Tchebycheff  polynomials  ,  if  we  know  the  value  of 
the  second  eigenvalue.   A-,,  for  then  all  the  eigenvalues 
except   A   lie  in  the  interval   (0,A-,),  and 

F(A)  =  Tjj(2|-  -  1) 

is,  as  is  well-known,  the  polynomial  of  degree  N  for 
which 

Max   |F(a)| 
0<A<A-,^ 

is  as  small  as  possible,  for  given  leading  coefficient. 
Tj^(z)   is  the  Tchebycheff  polynomial  defined  by 

The  zeros  of   T^(z)   are  clearly  at 

z  =  cos  ^1^  TT   ,     i  =  1,  2,  •••,  N  , 

and  therefore  we  take 

(16)  |i^  =  2^  (1  +  cos  ^1^  tt),  i  =  1,  2,  •••,  N. 

Since   A-,   is  not  known  in  advance,  we  use  an  approxi- 
mate value  of  it  which  is  improved  as  the  calculation 

see  Courant-Hilbert,  ibid.  Chapter  II,  section  9,  part  2. 
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proceeds,  by  a  calculation  of  \-.      and  the  eigenf unction 

X, (r)   that  proceeds  concurrently  with  the  calculation 

of   A   and  <l>(r).   The  calculation  of  X,  (r)   Is  similar 
o  1  * 

to  that  of  $(r),  except  that  it  is  necessarily  somewhat 
less  refined,  because  we  assume  no  knowledge  whatever  of  Xp. 
After  each  set  of  N  iterations  (15)  for  improvement  of  i>{r) , 
a  set  of  M  similar  iterations  is  performed  for  improvement 
of  X, (r).   These  are  of  the  same  form  as  (l^);  viz., 

i<p+l(r)  =  /  K^(r,s)kp(s)ds  -  v^kp(r) 

(17) 

W^^  ^  ^   ^a^^'^^VM-l^"^^"  -  Vp+M-l(^) 

except  that  the  constants   v   are  now  chosen  as 

V   =  A 
1    o 


(the  value  of   X   used  here  is  the  estimate  obtained 

^  o 

from  the  preceding  N  iterations  to  improve  ^(r)).   The 
corresponding  polynomial,  G(>v)  =  i'K-^^)  {'K-^^'^  "  •  {"K-^^)    is 
then  simply   (a-Aq)a^"'^.   It  vanishes  at   A  =  A^,  and  if 
one  is  lucky  in  the  choice  of   M  it  will  have  a  maximum 
near   A  =  X-,   and  be  quite  smaJ.l  for   A  =  A2J  etc.   (If 
one  does  not  wish  to  count  too  much  on  luck,  it  is  better 
to  have  M  too  large  than  too  small).   The  corresponding 
functions  k  (r)   approach  Xj^(r),  for  p— >oo,  as  ortho- 
gonality to  $(r)  is  assured  in  the  limit  because  G(>Vq)  =  0' 
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After  each  set  of  N  iterations  on  ^{r)      and  M 
Iterations  on   X-^(r),  the  value  of   a  is  adjusted  auto- 
matically in  such  a  direction  as  to  bring   A^   closer 
to   1.   To  do  this,  the  machine  remembers  the  last  two 
values  of  a  used,  say  a'   and  a",  and  the  correspond- 
ing approximations   A^  and  AjJ  to   X^.   It  then  takes 
as  the  next  trial  value  the  quantity 

a"(l-A^)+a'(A^-l) 

(18)      a=  ^TTTYT • 

o   o 

The  machine  then  prints  out  the  values  of   a  and  of  "K-^ 
obtained  so  far,  computes  the  new  kernel   K^(j^'S)  (in 
partially  integrated  form  as  described  below) ,  and  pro- 
ceeds with  the  next  iteration.   The  iterations  may  be 
stopped  (on  a  breakpoint)  whenever  desired:   the  values 
of  <t>{r)      and  X-,(r)  are  then  available  on  magnetic  tapes 
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IX.   The  Numerical  Procedure  for  Solving  the  Integral  Equation 

The  Integration  Is  performed  numerically  after  inte- 
gration by  parts  to  remove  the  singularity  of  the  kernel 
K  (r,s)   at   r  =  s.   Tnat  is,  we  set 


Ax  =  a/J,      (J  =  Integer); 


then.  In  order  to  evaluate  an  Integral  of  the  form 

a 


f(x)  =  /  K^(x,y)Q9(y)dy 


we  write 


J-1 


(19)         f,,,/2=^U) 


f   (j+l)Ax 

/       K  (x  ,,y)dy 
jAx      ^     1+1 


J  0+2" 


where  f.  -,  /p   is  an  abbreviation  for  f(  (1+1/2)  Ax) ,  and 
QP .  ,  /„   ^°^    ((j+1/2)Ax).   The  quantity  in  the  square 
bracket  is  a  coefficient  depending  on  i  and  J :  it  can 
be  computed  in  terras  of  the  function 

X 

g(x)  =  /  E(y)dy  , 
o 

a  table  of  which  has  been  computed  in  advance  by  a 
special  routine  and  is  available  on  tape. 

At  the  beginning  of  the  routine  there  is  an  oppor- 
tunity for  the  operator  to  type  in  values  of  J   (the 
number  of  net  points),  of  N  and  M   (the  numbers  of 
iterations  of  each  kind)  and  of  a   (the  radius  of  the 
sphere  in  mean  free  paths).   The  machine  then  generates 
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some  rather  crude  trial  functions  (slnusordal)   <I)(r)   and 
X-|(r)   and  then  goes  into  the  Iterative  routine. 

To  Illustrate:  with  J  =  kO,  N  =  M  =  3,  a  =  2.0, 
the  successive  iterations  gave  values  in  the  following 
table: 

a  A, 


-.^67 

-.500  .521 

-.4923  .522 

-.4925  .5223 

-.49250  .5222 

-.49249  .52222 

-.492482  .522210 

The  convergence  is  seen  to  be  quite  rapid;   (we  have  given 

only  those  digits  that  seem  relevant).   To  test  the  effect 

of  changing  the  interval  size  Ax,  a  second  run  was  made 

with  J  =  20  instead  of  40.   This  gave  a  similar  sequence 

of  numbers  appearing  to  tend  asymptotically  to  the  values 

a  =  -.49252   ,   Xj_  =  .52168 

It  was  on  the  strength  of  these  results  that  we  took 
a  =  -.4925  in  section  VI. 

The  functions  ^(r)   and  X-j^(r)   obtained  for  j  =  40 
are  shown  graphically  in  figure  7. 
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Figure  7 

In  another  series  of  calculations,  for  a  sphere  of 
^.0  mean  free  paths  radius,  i.e.  with  a  =  4.0,  the 
asymptotic  values  were 


J  =  20 
J  =  40 
J  =  60 


a  =  -.1551^ 
a  =  -.13295 
a  =  -.15289 


■K-^  =  .7557 
>^■^  =  .7564 
x-L  =  .7369 
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List  of  Machine  Routines 


A.   "El"  routine  computes  a  table  of  the  function 

X 

r 

o 


10  "six)    =  10"^  /  E(y)dy 


from  series  expansions  and  other  known  formulas, 
and  writes  a  table,  at  arguments: 

X  =  0  -  1.04  (.02), 
X  =  1  -  5.2  (.1)  , 
X  =  5  -  26   (.3)  , 

In  three  consecutive  blocks  on  a  magnetic  tape. 
The  degree  of  overlap  was  designed  to  facilitate 
quadratic  Interpolation.   There  Is  no  Input. 

B.  "Eep"  routine  prepares  the  preliminary  trial  values 
of  <I)(r)   and  X,  (r)  and  an  estimate  of  a  based 
on  slightly  modified  diffusion  theory.   Input:  a, J. 

C.  "Milne"  routine  computes  the  coefficients  in 
formula  (19)  for  the  numerical  Integration  from 
the  given  value  of  a  and  the  table  of  g(x), 

and  performs  the  iterations  (15)  and  (17)  described 
above.   For  the  first  two  sets  of  iterations   a 
is  taken  as  0.9a   and  l-loi   where  a   is  the 
preliminary  estimate  of  a  obtained  in  the  "Eep" 
routine.   After  that,  equation  (18)  is  always  used 


58- 


NYO-7696 

Physics 
to  obtain  the  new  a.   Input:   N  (M  Is  taken 
=  N  by  this  routine). 

D.  "Transport  eqn.  coef .  gen. "  prepares  In  advance 
and  writes  on  tape  two  sets  of  coefficients: 
one  for  Interpolating  among  the  tabulated 
values  (1)^  to  obtain  't>^{^-  .-,  /o    ■; )  ^^^  '^^^ 
for  the  double  interpolations  among  the  y}    . 
needed  in  computing  the  Integral  expression 
for  $(r)   by  equation  (9).   The  ability  to 
precompute  these  coefficients  and  read  them 
from  tape  as  needed  without  loss  of  computing 
time  proved  to  be  a  marked  improvement  of  the 
present  code  over  the  Los  Alamos  Maniac  code, 
which  had  been  prepared  Defore  the  magnetic 
drum  was  available  on  that  machine.   Input: 

I  (the  routine  takes  J  =  I) . 

E.  "Transport  eqn.  main"  (mod  2):  performs  the 
calculation  described  in  section  IV,  modified 
as  described  in  section  V.   After  each  cycle 

a  table  of  ^""^   and  values  of  various  integral 
quantities  are  printed  on  the  supervisory 
control  printer.   Various  breakpoint  options 
are  available;  e.g.  to  Increase  or  decrease 
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the  number  of  Iterations;  to  write  a  complete 
table  of  "i        (x,y)   at  any  desired  cycle  from 
the  working  tapes  onto  another  tape  where  It 
can  be  saved  for  subsequent  printing.   Three 
working  tapes  are  used:   In  any  cycle  one  is 
being  read,  one  Is  being  written  on,  and  one 
Is  Idle  (to  make  reruns  possible) :   The  three 
are  permuted  cyclically  after  each  complete 
cycle.   A  memory  dump  Is  made  every  cycle  for 
reruns.   Input:   a   (radius  of  sphere  In  mean 
free  paths) . 
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